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Abstract 

It is well-known that box filters can be efficiently computed using 
pre-integrations and local finite-differences [H |31 H] . By generalizing this 
idea and by combining it with a non-standard variant of the Central 
Limit Theorem, a constant-time or 0(1) algorithm was proposed in 
[3] that allowed one to perform space-variant filtering using Gaussian- 
like kernels . The algorithm was based on the observation that both 
isotropic and anisotropic Gaussians could be approximated using cer- 
tain bivariate splines called box splines. The attractive feature of the 
algorithm was that it allowed one to continuously control the shape 
and size (covariance) of the filter, and that it had a fixed computational 
cost per pixel, irrespective of the size of the filter. The algorithm, 
however, offered a limited control on the covariance and accuracy of the 
Gaussian approximation. In this work, we propose some improvements 
by appropriately modifying the algorithm in [3]. 

Keywords: Linear filtering, Gaussian approximation, Central limit 
theorem. Anisotropic Gaussian, Covariance, Box spline, Cartesian grid. 
Running sum, 0(1) algorithm. 

1 Introduction 

A space-variant filter is one where the shape and size of the filter is allowed 
to change from point-to-point within the image. Unlike the more standard 
convolution filter, the efficient computation of space-variant filters (diffu- 
sion filters) remains a challenging problem in computer vision and image 
processing [9]. It has been shown that efficient algorithms for space- variant 
filtering can be designed using spline kernels, particularly when the space- 
variance is in terms of the scale (or size) of the kernel. For instance, Heckbert 
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proposed an algorithm for adaptive image blurring using tensor-products 
of polynomial splines, where the image is filtered with kernels of various 
scales using repeated integration and finite-difference [4]. Based on similar 
principles, namely, the scaling properties of B-splines, Munoz-Barrutia et 
al. have developed an algorithm for fast computation of the continuous 
wavelet transform at different scales [10]. The idea was extended in [12] to 
perform space- variant filtering using Gaussian-like functions of arbitrary size 
. This was done by approximating the Gaussian using separable B-splines. 
The flip side of using a separable construction, however, was that it offered 
limited steerability and ellipticity. More recently, it was shown in [3j that 
this limitation could be fixed using certain non-separable splines called the 
box splines [7\. This, however, did not solve the problem completely as the 
associated filtering algorithm offered only a limited control on the accuracy 
of the Gaussian approximation. In this paper, we address these algorithmic 
limitations and provide some simple solutions. To introduce the problem and 
to fix the notations, we briefly recall the main results in While we keep 
the presentation as self-contained as possible, we have avoided reproducing 
all the technical details, and refer the readers to original paper for a more 
complete account. 



1.1 Background 

The key idea behind the fast 0(1) algorithm for Gaussian filtering in [3j is 
the use of a single global pre-integration, followed by local finite differences. 
At the heart of this is the following non-standard form of the Central Limit 
Theorenfl 

Theorem 1.1 (Gaussian approximation, [3]). For a given integer N > 2, 
let dk = {k- 1)tt/N forl<k<N. Then 

N^oo n ^^^^['^y J^i^^^^^k + ysm0k) \ = exp - y (x^ + y^)^ (1) 

This result tells us that the Gaussian can be approximated by the product 
of rescaled sincj^ functions that have been uniformly spread out over the 
half-circle. Note that each term in the product in ([T]) is obtained through a 
rotation and rescaling of sine (u) 1(f), the tensor-product between the sine 
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function and the constant function of unit height. The main idea behind 
([l]) is that sinc(t) behaves as 1 — + O(t^) close to the origin, so that 
[sinc(t/\/iV)]'^ ~ exp(— when N is large compared to t. 

So how exactly does this observation lead to an 0(1) algorithm for 
approximating Gaussian filters? By 0(1), we mean that the cost of filtering 
is independent of the size of the filter. It is clear that the right side of ([T]) is 
simply the Fourier transform of the Gaussian 
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On the other hand, the product in ([T]) can be related to the box distribution 
BoXa(3;)5(y), where 



Box„(t) 



1/a for -a/2 < t < a/2, 
otherwise, 



and 5{y) is the delta function. Indeed, using the rotation-invariance and the 
multiplication-convolution property of the Fourier transform, we see that the 
left side of ([T]) is the Fourier transform of 

Boxai(xcos0i + ysm6i)5{xs\ii6i — ycosOi) * • • • 

BoXaj^(a;cos^Ar + ?/sin07v)'^(3;sin6'7v — y cos^at), (2) 



where every a^ equals a^JlA^jN . This gives us the dual interpretation of ([!]) 
in the space domain, namely that g{x^ y) can be approximated using the 
convolution of uniformly rotated box distributions. This idea is illustrated 
in Figure [!} The function in ^ turns out to be a compactly-supported 
piecewise polynomial of degree < N — 2, popularly called a box spline. Note 
that we could as well replace the sine function in ([T]) by any other "bump" 
function that looks like an inverted parabola around the origin, e.g. the 
cosine function as used in [§]. However, not every such function would be 
the Fourier transform of a simple, compactly-supported function such as the 
box function. 

Applied to a continuously-defined function f{x,y), this means that we can 
approximate the Gaussian filtering {f*g)(x, y) by successive convolutions with 
a fixed number of rotated box distributions. Now, note that the convolution 
of f{x, y) with BoXa{x cos Oj^ + y sin 9k)S{x cos O^ — y sin 0^) amounts to doing 
one-dimensional convolutions of BoXa(t) with profiles of f{x, y) sampled along 
Lq = {(x, y) : x cos Ok — ysinOj. = 0} and lines parallel to Lq. For example, 
the convolution of f{x,y) with Boxa(x)(5(y) is given by one-dimensional 
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Figure 1 : Demonstration of the Central Limit Theorem used to approximate 
the isotropic Gaussian, (a) Convolution of four box distributions of equal 
width and uniformly distributed over [0,7r), (b) Convolution of eight box 
distributions of equal width and uniformly distributed over [0,7r), (c) Target 
Gaussian. The widths of the boxes are so adjusted that (a) and (b) have the 
same covariance as (c). Note that the convergence happens quite rapidly 
with the increase in the number of boxes. The maximum pointwise error 
between (b) and (c) is already within 1% of the peak value. 

convolutions of Boxa(t) with the profiles of f(x,y) along the x-coordinate 
and lines parallel to it. The idea in [3J was to decompose the convolution 
with a box distribution into two steps: 

• Pre-integration of the scan profiles along Lq^ (1 < A; < N) and lines 
parallet to it. This resulted in the so-called running sum. 

• Application a finite-difference mesh to the running sum at each pixel, 
where the mesh parameters were determined by the scales ai, . . . , oat. 

By doing so, the cost of the convolutions with the box distributions was 
reduced to the cost of (1) and (2), which is clearly 0(1) per pixel irrespective 
of the values of oi a^r . 

This idea was then extended to perform the more challenging space- 
variant filtering. This was based on the observation that by adjusting the 
scale of the box spline along each direction 9i, ... ,6k, one could control 
the anisotropy of the box spline. From the algorithmic perspective, this 
meant that we could change the filter at every pixel simply by changing the 
finite-difference mesh. This resulted in a fast algorithm for space-variant (or, 
non-convolution) filtering that had the same 0(1) complexity. 
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In particular, for the case = 4, corresponding to 9i = 0, 02 = n / A, 9^ = 
7r/2, ^4 = 37r/4 in ([2]), a simple algorithm was developed that allowed one to 
control the covariance simply by adjusting the scales oi, . . . , 04 in (§. We 
continue to use /3a (x, y) to denote this four- directional box spline, where we 
call a = (ai, a2, as, 04) the scale-vector. The algorithm for the four-directional 
box spline proceeded in two steps. In the first step, the running sum was 
obtained by pre-integrating teh image along the directions 0,7r/4, 7r/2, and 
37r/4. This was done using fast recursions. To keep the running sum on 
the Cartesian grid, a step size of unity was used along the horizontal and 
vertical directions, while a step of was used along the diagonals. In 
the second step, a 16-tap finite-difference mest was applied at point of the 
running sum, where the taps of the mesh was determined by the scale-vector. 
This produced the final filtered image. We refer the readers to Algorithm 1 
in [3] for further details. Henceforth, we will refer to this as the BoxFilter 
algorithm. More recently, similar ideas have also been applied in the field 
of computer graphics [6]. By adapting the same algorithm, an efficient 
technique for detecting blobs in cell images was also proposed in 

The anisotropy of /3a(x,y) was specified using three parameters, its size 
(sa), elongation (pa), and orientation (^a)- These were defined using the 
eigen decomposition of the covariance matrix of f3a{x,y), 

f f x'^Pa{x,y) dxdy f f xyf3a{x,y) dxdy \ 
J J xy(3a{x,y) dxdy J f y'^f3a{x,y) dxdy J' 

The size was defined as the sum of the eigenvalues, the elongation as the 
ratio of the larger to the smaller eigenvalue, and the orientation as the angle 
(between and vr) of the top eigenvector. In fact, the covariance of /3a{x, y) 
was explicitly computed to be 



1 / 2af + a| + a\ 



This resulted in the following formulas in terms of the scale-vector: 



(4) 



(5) 



1 2 E "I + , n 4-al + VD 

where D = (af — 03)^ + (03 — a\)'^. 

The scale vector was to control the anisotropy of the filter at every 
pixel, and this was done, in effect, by controlling the tap positions of the 
finite-difference. Note that Ca = (a^/6)I when all the four scales are equal 
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to a. As a result, the best approximation of the isotropic Gaussian with 
variance cr^ was obtained by seeting each scale to ^/Ea. As for the anisotropic 
setting, a simple algorithm for uniquely determining the scales for a given 
size, elongation, and orientation (equivalently, for a given covariance) was 
provided in [3]. The algorithm took the covariance C as input and returned 
the box spline that had the minimum kurtosis among all box splines with 
covariance C. 



1.2 Problem description 

There are, however, two bottlenecks with the above approach. The first is 
that one cannot control the accuracy of the Gaussian approximation in case 
of the space-variant filtering. We note that this problem can be addressed 
for the more simple setting of convolution filtering simply using higher-order 
box splines (see Figure [l]). In effect, this amounts to repeated filtering of the 
image with the four- directional box spline However, this cannot be done in 
the space- variant setting, where the kernel changes from point-to-point. The 
only possibility is to adapt the algorithm for the higher-order box splines. 
This, however, turns out to be of limited practical use since the size of the 
finite-difference mesh grows as 2^ with the order (two points per box 
distribution). 

The other problem is the control on the elongation of the anisotropic 
box splines. This problem is intimately related to the geometry of the 
Cartesian grid Z'^. Note that the axes Lq of the box distribution in the 
continuous setting corresponds to the discrete points Lq = {{m,n) G : 
m cos 9k — nsinOk = 0} on the grid. As a result, one requires tan^fc to be a 
rational number to avoid interpolating the off-grid samples. It is, however, 
easily seen that, for A > 4, one cannot find lines Lq-^, . . . ,Lq^ with the 
requirement that 9i, . . . ,9n are uniformly distributed. This restricted us to 
the four-directional box spline (3aix,y)- This restriction on the number of 
directions resulted in a ceiling on the maximum achievable elongation. In 
particular, while Paix,y) could arbitrarily elongated in the neighborhood of 
its four axes, it was difficult to achieve large elongations in the neighborhood 
of the mid-axes, i.e., along vr/S, Svr/S, 57r/8, and Ttt/S. It was found that, 
for every orientation < (j) < n, there is a bound e{(j)) on the elongation 
achievable by /3a (x, y) given by 



e(0)=^ + ^ + f±g, (6) 
where t = \ tanc/) — cot(/)|/2. This meant that pa could be made at most 
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as large as e{(j)) under the constraint that 9^ = 4>, where and 9^ are as 
defined in ([s]). The smahest value of e(0) as (p varied over the half-circle was 
found to be roughly 5.8, and this was reached along the mid- axes where one 
had the least control. 

2 Improved accuracy and elongation 

We now address the above mentioned problems for the four-directional box- 
spline. Our interest is mainly in the space-variant setting, where the filter 
varies from point-to-point in the image. 

2.1 Improved accuracy 

As remarked earlier, the Gaussian approximation can be improved by con- 
volving the box spline with itself. All that needs to be guaranteed is that the 
covariance remains equal to that of the target Gaussian after the convolution. 
In fact, when the functions involved have some structure, one can exactly 
predict how the covariance changes with smoothing. 

Proposition 2.1 (Covariance of convolutions). Let f{x,y) and g(x,y) be 
two functions which are symmetric around the origin, and which have unit 
mass. Then the covariance of (/ * g){x, y) is the sum of the covariances of 
f{x,y) andg{x,y). 

Note that f3a{x, y) automatically satisfies the conditions in the proposition. 
Our next observation is that we can get a higher-order space- variant filter 
from a lower-order space- variant filter by applying a global convolution. We 
note that the space- variant filtering of /(x,y) can generally be written as 



where h{x, y; n, v) is the space- variant kernel that is indexed by the spatial 
coordinates {x,y). 

Proposition 2.2 (Filter decomposition). Let hi{x, y; u, v) be a space-variant 
kernel, and let h2{x,y) be a convolution kernel. Then the convolution of 
f{x, y) with h2{x, y), followed by the space-variant filtering with hi{x, y; u, v), 
can be expressed as a single space-variant filtering: 
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where 



h{x,y;u,v) 



oo roo 



hi{x, y; u , v')h2{u — u ,v — v') du'dv. 



oo J — oo 



Propositions 2.1 and 2.2 together suggest a simple means of improving 
the Gaussian approximation for space-variant filtering. Suppose that C is 
the covariance of the target Gaussian. We split this into two parts, namely, 



C = cj^I + AC. 



(7) 



Note that C and AC have the same orientations, that is, cr^I possibly alters 
the elongation, but not the orientation. This is confirmed by the fact that C 
and AC have the same set of eigenvectors. 



Algorithm 1 Space-variant 0(1) Gaussian filtering (better accuracy) 
Input: Image f{m,n), covariance map C(m, n). 

1. Set (T^ to be half the bound in ([9]). 

2. Set ai = • • • = 04 = \/6cT. 

3. Pass ai, . . . , 04 and /(m, n) to the BoxFilter algorithm, which return 
g{m,n). 

4. At each pixel, compute a(m,n) from C{m,n) — cj^I. 

5. Pass a(m, n) and g{m, n) to BoxFilter algorithm, which returns /(x, y). 
Return: Filtered image f{x,y). 



The idea is we first convolve the image f{x,y) with an isotropic (3{x,y) 
of covariance cr^I (each is ^/Qa). This is done using the 0(1) algorithm 
described earlier. We then calculate the residual AC at each point in the 
image. This is used to fix a of the anisotropic I3g,{x, y). The scale assignments 
are then used for the space- variant filtering, again using the 0(1) algorithm. 
The main steps are given in Algorithm [T] It is clear that the overall algorithm 
requires 0(1) operations per pixel. 

There is a technical point that must be addressed at this point. Namely, 
cr must be so choosen that AC is again a valid covariance matrix, that is, AC 
be positive definite. Now, note that AC can be written as Q^(A — cj^I)Q, 
where Q-^AQ is the eigen decomposition of C. It follows that a necesaary 
and sufficient condition for AC to be a valid covariance matrix is that o"^ 
must be smaller than the minimum eigenvalue of C, that is, 

cJ^<l (Cii + C22 - ^(Cii-C22)2 + 4Cf2] . (8) 
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Figure 2: This show the steps involved in improving the accuracy of the 
Gaussian approximation. The isotropic ones are shown on the top panel, 
while the bottom panel shows the anisotropic ones. In either case, (a) is the 
four-directional box spline with covariance and (b) is the four-directional 
box spline with covariance C — cr^I, where C is the covariance of the target 
Gaussian. The resulting box spline approximation is shown in (c), which 
indeed looks more Gaussian- like than the ones in (b) and (c). See Figure [3] 
for a comparison of the approximations. 
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The other consideration comes from the bound in Q (this concerns only 
the anisotropic case), which requires that the the ratio of the eigenvalues of 
AC must be with e{(j)), where (p is the orientation of the target Gaussian. A 
simple calculation shows that this is equivalent to 

< \ (Cn + C22 - (^If^) V(Cii-C22)2 + 4C?2) . (9) 

Since e{(j)) > 1, this bound is clearly smaller than the one in ([s]), so it suffices 
to guarantee that ([9]) holds. 

We have empirically observed that, while it is the order of the box splines 
that dictates the final smoothness, it is the size of a that controls the level of 
approximation. In practice, we have observed that the best approximations 
are obtained when o"^ is roughly 50% of the bound in Q . We will discuss 
this in detail in the sequel. 

2.2 Improved elongation 

We now address the elongation problem. The elongations achievable by 
the box spline is minimum exactly mid-way between the axes of the box 
distributions, namely along 22.5°, 67.5°, 112.5°, and 157.5°. An intuitive 
solution to this is to place box distributions along these angles. As mentioned 
earlier, the problem is that the tangents of these angles are not rational. The 
best we can do in this situation is to approximate the tangents by rational 
numbers. In other words, for an given angle 0, the problem reduces to one of 
finding integers p and q such that p/q ^ tan 6*. Of course, we would also like 
to keep p and q as small as possible. Also, for every 9 with tan 6* = p/q, we 
balance it with the orthogonal angle 6 + 90°. This is easily done by setting 
the tangent of the orthogonal angle to —q/p- For example, if we set p = 1 
and q = 2, then we get two orthogonal directions at roughly 26.6° and 116.6°. 
The other two directions are obtained by setting p = 2 and q = I which gives 
the orthogonal directions at roughly 63.4° and 153.4°. Though these are 
not distributed uniformly over the half-circle, they are off the mark by only 
about 4°. To the best of our knowledge, this is the best we can do keeping p 
and q as small as possible. 

This gives us eight directions on the Cartesian grid that are "almost" 
uniformly distributed over the half-circle and whose tangents are rational. 
The idea then is to place box distributions along all or some of these directions. 
The higher the number of directions, the better is the Gaussian approximation. 
Unfortunately, as mentioned in the introduction, there are practical difficulties 
in simultaneously using the eight directions. Keeping in mind that we 
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already have a solution to improve the accuracy, we instead propose the 
following trick. The idea is to use a pair of four-directional box splines, 
/3a(x,y) with angles G = {0°,45°,90°, 135°}, and P'^{x,y) with angles 6' = 
{26.6°, 63.4°, 116.6°, 153.4°}. We next divide the half-circle into two disjoint 
sectors S and S' , where S consists the neighborhood of the former angles, 
and S' consists of the neighborhood of the latter angles. To approximate a 
Gaussian with orientation cp, we use I3a{x,y) if (p is in S, and (3'g^{x,y) if </) 
is in S' . This allows us to achieve better elongations than what /3a (2;, y) or 
/3a(a;, y) could achieve by themselves. Of course, we will now require two sets 
of running sums, one corresponding to each of the box splines. The overall 
algorithm for space- variant filtering is summarized in Algorithm [2j We note 
that, as proposed in Section 2J_, we could smooth the four-directional box 
splines with a isotropic box spline, if higher accuracy is desired. 



Algorithm 2 Space-variant 0(1) Gaussian filtering (improved elongation) 

1. Input: Image f{m,n) and covariance map C(m,n) on Z^. 

2. Perform running sums on f{m,n) along Q to get g{m,n). 

3. Perform running sums on f (171,12) along Q' to get g'{m,n). 

3. At each (m, n) do 

(a) Compute the orientation (j) from C(m, n). 

(b) if (/> e 5, then 

(bl) Compute a for /3a(a;,y), and the finite-difference parameters. 
(b2) Apply finite-difference on g{m,n) to get f (171,71). 
else (0 G S') 

(bl) Compute a for (3'^{x,y), and the finite-difference parameters. 
(b2) Apply finite-difference on g'{m,n) to get f{m,n). 

4. Return: Filtered image f{m,n). 



2.3 Filtering algorithm for f3'^{x,y) 

We note that the algorithm for space-variant filtering using pi^{x, y) is identi- 
cal to the one for f3a{x,y). For completeness, we describe this in Algorithm 
[3] Here xq, . . . , X15 and yo, . . . , yi5 are the coordinates of the finite-difference 
mesh, and wq, ■ ■ ■ , W15 are the corresponding weights. The exact values are 
given in table 1. The shifts ri and T2 along the image coordinates are given 
by 

Ti = — ^(2ai+a2 — 03 — 204) and r2 = — ^(ai+2a2+2a3+a4— 6\/5). (10) 
2v5 2v5 
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Table 1: The positions and weights of the finite-difference mesh used in 
Algorithm [sj Here a = (01,02,03,04) are the scales of the mesh, and 
w = 1/(01020304) and a'- = ai/y/5. 
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(0,0) 


+w 
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(-204,04) 


—w 
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(2o'i,o;) 


—w 
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- 204, a[ + 04) 


+w 
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(o'2,2o'2) 


—w 


10 
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- 20^1,20; + 0^1) 


+w 


3 


(2oi 


+ a'2, a'l + 202) 


+w 


11 


(2a'i 


+ o'2 


-2o;t,o; + 2o'2 + o'4) 


—w 


4 




(-o'3,2o'3) 


—w 


12 




(-4 


-2a'4,2o'3 + o'4) 


+w 


5 


(2o; 


-a'3,a; + 2a'3) 


+w 
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(2o'i 


-o'3 


-2o;t,o; + 2o'3 + o'4) 


—w 


6 


(4- 


-a'3,2a'2 + 2a'3) 


+w 


14 


(4- 


-«3 - 


2o^,2o2 + 2o'3 + 0^) 


—w 


7 


(2a'i + 02 


- 03,0; + 202 + 203) 


—w 


15 


(2o; + 02 


-03- 


- 204, a[ + 202 + 203 + 04) 


+w 



Note that we have set the scales of the running to \/5 along each of 
the four directions. This ensures that no off-grid sample, and hence no 
interpolation, is required to compute the running sums. The sample F{x, y) 
used in the finite-difference is given by 

F{x,y)= ^ g4.{m,n)P[,{x-m,y-n) (11) 

\m—x\<3, |n— 3/|<3 

where b = {V^, \/5, \/5 , -y/S) . The box spline (3[^{x,y) is supported within 
the square [—3,3]^, and this is why we have a finite sum in (11). This can 



be seen as an interpolation step, where we interpolate the running-sum using 
the samples of the box spline. 

2.4 Analytical formulas for /3^(x,y) 



The exact formula for the box spline /3(,(x,y) in (11) can be computed 
analytically. To do so, we note that this box spline can be expressed as the 
convolution of two rescaled and rotated box functions. In particular, (3^{u,v) 
is given by the area of overlap between a box function of width \/5 and height 
l/\/5 that is centered at the origin and rotated by (f>i, and a second box 
function of same width and height that is centered at (u, v) and rotated by 
02, where tan(/)i = 1/2 and taii(p2 = 2. The overlapping region is a polygon, 
whose area can be determined by inspection. This can also be found using 
symbolic computation. For example, this can be done in Mathematica by 
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setting 

f[x,y] := (l/\/5) UnitBox[(ucos(/)i + vsin0i)/\/5] 
g[x,y] := (l/\/5) UnitBox[(ucos(/)2 + vsin02)/\/5] 

and then invoking the command Convolve[ f[u, v], g[u, v], {u, v}, {x,y}]. 
The function returned is a piecewise polynomial of second degree that is 
compactly supported on a octagonal domain. The equation for each patch is 
stored in a look-up table. 



Algorithm 3 Space-variant 0(1) filtering using f3'i^{x, y) 

1. Input: Image f{m,n) and scale map a(m, n) on Z^. 

2. Use recursion to compute the following: 

(a) gi{m,n) = \/5f{m,n) + gi{m -l,n- 2). 

(b) g2im,n) = \fhgi{m,n) + g2{m - 2,n - 1). 

(c) g3im,n) = ^/bg2{m,n) + g3{m + l,n- 2). 

(d) g4{m,n) = V5g3{m,n) + g4{m + 2,n- 1). 

3. At each (m, n) do 

(a) Set up Wi,Xi,yi{i = 0, 1, ... , 15) using table [T| and ri, T2 using (10). 

(b) Compute Fi = F{m + ti — Xi,n + T2 — yi) for i = 0, 1, . . . , 15 using 
©• 

(c) Set f{m,n) = wqEq + ■■■_ + wi^Fi^. 

4. Return: Filtered image f{m,n). 



Similar to f3aix,y), we can control the covariance (size, elongation, and 
orientation) of (3i^{x,y) using its scale vector. We computed the covariance 
of P'g,{x, y) to be 
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4af + al + al + Aal 2(af + a\ 



2{af + a: 



a|) af + 4a^ + 4a| + a|) 



We defined the size, elongation and orientation as in ([s]). The final expressions 
are respectively 



r 

60' 



where p = 3(af 
al + al + al). 



and tan 



Y^p2 _|_ g2 ' 



-p + Y^p2 _|_ ^2 



al), and r = 5(af -|- 



We note that ai, . . . , 04 cannot be uniquely determined from the above 
three constraints. As proposed in our earlier paper, we select the box spline 
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that has the minimum kurtosis of ah box sphnes with a given covariance. 
This turns out to be a one-dimensional optimization problem and can be 
solved efficiently using a root-finding algorithm. The kurtosis of P'^{x,y) 
using the definition in [13j turns out to be 

K / 4a| + a| + a| + 4aj 2(af + a| - a| - a|) \ 
5 V 2(af + aj-4- aj) af + 4a| + 4a| + a^) J' 

where k is the kurtosis of the box function Boxi(t). To get a number from 
this matrix, we simply take its Frobenius norm, which is then optimized 
to get the scale vector. The steps of the derivation are identical to the one 
described in [3], and therefore we skip the details. 

3 Experiments 

We now present some results that confirm the improvements suggested in the 
previous section. First, we try to see the improvement in the accuracy of the 
Gaussian approximation. To do so, we consider an isotropic Gaussian with 
C = I, and an anisotropic Gaussian with size 1, elongation 3, and orientation 
30°. In either case, we split the covariance as C = cj^I -|- AC. The steps 
in the construction are shown in Figure Q, and the final box splines are 
compared in Figure ([3]). In either case, we have used the four-directional 
box spline to approximate the Gaussian with covariance a^I, where the four 
scales were set to y/6a. To approximate the Gaussian with covariance AC, 
we used a four-directional box spline, but now with unequal scales. The 
four scales were computed by minimizing the kurtosis of the box spline as 
explained earlier. It is seen from the figures that the final box splines are 
indeed very good approximations of the target Gaussian. In particular, the 
maximum pointwise error was within 1% of the peak value for the isotropic 
one, while it was within 2% for the anisotropic ones. Note in figure ^ how 
the convolution helps in "rounding" the sharp corners of the four-directional 
box spline. 

To get a better understanding of the improvement in the accuracy, we 
compare the approximation error between the target Gaussian and the two 
box splines. We do this for various sizes, elongations, and orientations of the 
target Gaussian g(x,y). We use the ratio ||/ — 5112/115112 as the normalized 
error between the approximating box spline f{x,y) and g{x,y), where ||/||2 
is the L2 norm of a function. We computed the error from the analytical 
formulas of the Gaussian and the box splines using numerical integration. The 
results are reported in table ^ In the first two instances (isotropic setting). 
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(a) (6) 



0\0\o 



(a) (6) (c) 

Figure 3: In this figure, we compare the target Gaussian (c) with the box 
sphne approximation (a) proposed in [3j and the one (c) proposed in this 
work. The isotropic Gaussian has covariance I, while the anisotropic one has 
elongation 3 and orientation 30°. The maximum pointwise error is within 
1% of the peak value for the former and within 2% for the later. A visual 
comparison confirms the improvement in the accuracy - it is indeed hard to 
distinguish (b) from (c). 



Table 2: Improvement in Gaussian approximation using Algorithm [TJ The 
table shows the normalized errors for the BoxFilter in |3j (second row) 
and the filter proposed in the paper (third row). This is done for certain 
representative size (s), elongation (p), and orientation (6) of the target 
Gaussian. 





(1,1,0) 


(5,1,0) 


(1,4,0) 


(5,4,0) 


(5,3,^/8) 


(5,8,7r/3) 


(5,5,7r/2) 


Error (%) in 


m 


10.8 


10.8 


19.1 


18.7 


23.9 


19.2 


17.2 


Present error 


(%) 


4.9 


4.9 


15.3 


14.6 


20.8 


15.8 


12.6 


Improvement 


(%) 


54.8 


54.8 


21.6 


18.32 


12.7 


17.7 


26.7 
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note that there is almost a 40% improvement in the approximation (this is 
the case for any size). The improvement is less dramatic, but nonetheless 
appreciable, for the remaining five instances (anisotropic setting). We found 
that this is between 10% — 30% for the anisotropic box splines in general. 
The improvement is on the higher side when the orientation is close to one 
of the four axes of the box spline. 

We now discuss the choice of ex in ([T]) used in the above experiments. 
As noted earlier, we have observed that the best results are obtained when 
a for the isotropic box spline is roug hly 50% of the bound in The 
variation in the approximation accuracy is shown in table [3] for a particular 
Gaussian. Exhaustive experiments across various settings have shown that 
this happens consistently across different size, elongation, and orientation. 
The explanation for this is as that very little smoothing takes place when 
a is small, which explains the marginal improvement in accuracy. On the 
other hand, large values of a result in over-smoothing, where the isotropic 
box spline dominates the overall anisotropy - the contribution from the 
anisotropic part is less significant. Following these observations, we set cr^ to 
be 50% of the bound in 



Table 3: Effect of cr in Q on the Gaussian approximation. The table shows 
the difference between (a) the approximation error between a fixed Gaussian 
(s = 5,/9 = 3,^ = 7r/4) and the BoxFilter, and (b) the error between the 
same Gaussian and the present filter. The best results are obtained when cr^ 
is roughly 50% of the bound in ([9|). 



% of the bound in ( 


9 


) 


10 


20 


30 


40 


50 60 


70 80 


Improvement (%) 


3.5 


12.1 


22.3 


29.9 


33.3 33.13 


30.9 28.8 



We now study the improvement in elongation obtained using Algorithm [2j 
While there is an explicit formula for the bound of /3a(x, y), we computed the 
bound of I3'^{x,y) empirically. Table |4] gives a comparison of the maximum 
elongations achievable using I3a{x,y) and j3'g^{x,y) as against that achievable 
using only /3a(a;, y)- In the table, we give the bounds on the elongation for 
orientations in the sector [0°,45°], spanned between two neighboring axes 
of j3a{x,y). While large elongations are achievable using Pa{x,y) in the the 
caps [0°, 5°] and [40, 45°], the bound falls off quickly outside these caps. In 
particular, the bound falls to a minimum of 5.8 at 22.5°. The second box 
spline f3'g^{x,y), with its axis along 26.6°, helps us improve the elongation 
over [5°, 40°]. For example, the bound at 22.5° has now increased to 10.8, 
and remains fairly high in the cap [20°, 40°]. However, a minimum of 8.2 is 
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reached at 13.3°, exactly midway between the axes of /3a{x,y) and j3'^{x,y). 
This is nevertheless better than the previous low of 5.8. We have evidence 
that the use of a third box spline with axes close to 13.3° can improve 
the minimum bound to 15.2. We note that the use of multiple box splines 
increases the run time of the BoxFilter by only 10% — 15%. The added 
computations are required only from the running sums, which is the fastest 
part of the algorithm. The finite-difference part of the algorithm remains the 
same. For every orientation, one simply has to choose the correct running 
sum and the associated finite-difference. 

Table 4: Comparison of the elongations achievable in [3] and using Algorithm 
[2| The table shows the bounds on the elongation at different orientations 
using a single box spline filter (second row) and a pair of box spline (third 
row). The new bound is the maximum of the bounds of the two box splines. 



Orientation 


0° 


5° 


13.3° 


20° 


22.5° 


25° 


26.6° 


30° 


40° 


45° 


Previous bound 


oo 


13.6 


6.8 


5.9 


5.8 


5.9 


6.0 


6.1 


13.6 


oo 


New bound 


oo 


13.6 


8.2 


9.5 


10.8 


30.1 


oo 


28.8 


13.6 


oo 



We have implemented a basic version of Algorithm [T] as an Image J plugin. 
The Java code can be found at www.math.princeton.edu/~kchaudhu/code.html 
Some of the results in this paper can be reproduced using this plugin. The 
typical run time on a 256 x 256 image (Intel quad-core 2.83 GHz processor) 
was around 300 milliseconds. As expected, the run time was approximately 
the same for filters of different shape and size. On the other hand, we have 
tested the formulas for the box spline /3a (x, y) by implementing Algorithms [2] 
in Matlab. The typical run time was around 600 milliseconds on a 256 x 256 
image. The results obtained on a impulse image are shown in Figure [4j In 
the future, we plan to improve the run time by implementing the algorithms 
in Java. 

4 Discussion 

In this work, we showed that the performance of the box spline filters in 
[3] can be improved, while preserving the 0(1) complexity of the orignal 
algorithm. In particular, we showed how the accuracy of the Gaussian 
approximation can be improved by pre-filtering the image with an isotropic 
Gaussian of appropriate size. This idea is in fact quite generic, and can be 
used to improve the performance of any algorithm for space-variant filtering. 
We also showed that better elongations can be achieved using a pair of 
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(a) (b) (c) 

Figure 4: This shows the Gaussian-Uke blobs obtained by applying Algorithm 
[3] to an impulse image, (a) is the isotropic response obtained using equal 
scales, (b) and (c) are the anisotropic responses obtained using unequal 
scales. The latter are oriented along 45° and 135°, and have elongation 3. 

four-directional box splines instead of one. The present work suggests that 
both the accuracy and the elongation can be dramatically improved, at a 
small extra cost, by using a bank of four-directional box splines (ideally, 
three or four) whose axes are evenly distributed over the half-circle. We will 
investigate this possibility in the future. 

5 Acknowledgment 

The authors would like to thank Professor Michael Unser and Dr. Arrate 
Munoz-Barrutia for valuable discussions. This work was supported by the 
Swiss National Science Foundation under grant PBELP2-135867. 

References 

[1] P. Viola and M. Jones, "Rapid object detection using a boosted cascade 
of simple features," IEEE Conference on Computer Vision and Pattern 
Recognition, pp. 511-518, 2001. 

[2] F.C. Crow, "Summed-area tables for texture mapping," ACM Siggraph, 
vol. 18, pp. 207-212, 1984. 

[3] K.N. Chaudhury, A. Munoz-Barrutia, M. Unser, "Fast space-variant ellip- 
tical filtering using box splines," IEEE Transactions on Image Processing, 
vol. 19, no. 9, pp. 2290-2306, 2010. 



18 



[4] P. S. Heckbert, "Filtering by repeated integration," International Confer- 
ence on Computer Graphics and Interactive Techniques, vol. 20, no. 4, 
pp. 315-321, 1986. 

[5] K.N. Chaudhury, "A radial version of the Central Limit Theorem," 
|arXiv:1109.2227i yl [cs.IT]. 

[6] G. Wetzstein, D. Lanman, W. Heidrich, and R. Raskar, "Layered 3D: 
tomographic image synthesis for attenuation-based light field and high 
dynamic range displays,", ACM Transactions on Graphics, vol. 30, no. 4, 
pp. 1-12, 2011. 

[7] C. de Boor, K. Hollig, and S. Riemenschneider, Box Splines, Springer- 
Verlag, 1993. 

[8] K.N. Chaudhury, "Constant-time filtering using shiftable kernels," IEEE 
Signal Processing Letters, 2011. 

[9] J. Weickert, "Theoretical foundations of anisotropic diffusion in image 
processing," Theoretical foundations of computer vision, vol. 11, pp. 
221-236, 1994. 

[10] A. Munoz-Barrutia, R. Ertle, and M. Unser, "Continuous wavelet trans- 
form with arbitrary scales and 0{N) complexity," Signal Processing, vol. 
82, no. 5, pp. 749-757, 2002. 

[11] K.N. Chaudhury, Zs. Puspoki, A. Munoz-Barrutia, D. Sage, and M. 
Unser, "Fast detection of cells using a continuously scalable Mexican- hat- 
like template," IEEE International Symposium on Biomedical Imaging, 
pp. 1277-1280, 2010. 

[12] A. Munoz-Barrutia, X. Artaechevarria, and C. Ortiz-de Solorzano, "Spa- 
tially variant convolution with scaled B-splines," IEEE Transactions on 
Image Processing, vol. 19, no. 1, pp. 11-24, 2010. 

[13] A. Tkacenko and P.P. Vaidyanathan, "Generalized kurtosis and applica- 
tions in blind equalization of MIMO channels," Conference Record of the 
Thirty-Fifth Asilomar Conference on Signals, Systems and Computers 1 , 
742 - 746, 2001. 



19 



